Figure 1 - This figures shows…
library(png)
library(maptools)
Checking rgeos availability: TRUE
library(raster)
library(gam)
Loading required package: splines
Loading required package: foreach
foreach: simple, scalable parallel programming from Revolution Analytics
Use Revolution R for scalability, fault tolerance and more.
http://www.revolutionanalytics.com
Loaded gam 1.14
setwd("~/Desktop/Botero postdoc 2016/Human density and the origins of agriculture/")
par(mar=c(0,0,0,20))
d <- readPNG("Larson_dates.png")
plot(seq(0,18, length.out = 19), seq(0,36, length.out = 19), type="n",ylim=c(0,36),xlim=c(0, 18), xaxt="n")
rasterImage(d, 0,0,18,36, interpolate=TRUE, col=d)
Start_of_early_window <- 16-12
End_of_early_window_start_of_late_window <- 8.2
End_of_late_window <- 17-4.2
polygon(x=c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window), y=c(0, 34, 34, 0), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x=c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window), y=c(0, 34, 34, 0), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
These dates are provided in the supplimentary information for the Larson (2014) paper. I’ve copied those values into a .csv table provided here.
domestication_times <- read.csv("Domestication timing larson 2014.csv")
dim(domestication_times)
[1] 77 8
| Region | Species | Start.Exploitation | Finish.Exploitation | Start.predomestication | Finish.predomestication | Start.Domestication | Finish.Domestication |
|---|---|---|---|---|---|---|---|
| Southwest asia | Wheat | 12.00 | 11.25 | 11.25 | 11.00 | 11.00 | 9.00 |
| Southwest asia | Barley | 12.00 | 11.25 | 11.25 | 10.50 | 10.50 | 9.00 |
| Southwest asia | Lentil | 12.00 | 11.00 | 11.00 | 10.50 | 10.50 | 9.00 |
| Southwest asia | Pea | 11.50 | 11.00 | 11.00 | 10.00 | 10.00 | 8.50 |
| Southwest asia | Chickpea | 11.00 | 10.50 | 10.50 | 10.25 | 10.25 | 8.25 |
| Southwest asia | Broadbean | NA | NA | NA | NA | 10.50 | NA |
| Southwest asia | Flax | 12.00 | 9.50 | NA | NA | 9.50 | NA |
| Southwest asia | Olive | 10.00 | 6.00 | NA | NA | 6.00 | NA |
| Southwest asia | Sheep | 12.00 | 10.50 | 10.50 | 9.75 | 9.75 | 8.00 |
| Southwest asia | Goat | 12.00 | 10.50 | 10.50 | 9.75 | 9.75 | 8.00 |
| Southwest asia | Pig | 12.00 | 11.50 | 11.50 | 9.75 | 10.25 | 9.00 |
| Southwest asia | Cattle, taurine | 11.50 | 10.50 | 10.50 | 10.25 | 10.25 | 8.00 |
| Southwest asia | Cat | NA | NA | 10.50 | 4.00 | 4.00 | NA |
| South Asia | Tree cotton | 8.50 | 4.50 | NA | NA | 4.50 | NA |
| South Asia | Rice | 8.00 | 5.00 | 5.00 | 4.00 | 4.00 | 2.50 |
| South Asia | Little millet | NA | NA | NA | NA | 4.50 | NA |
| South Asia | Browntop millet | NA | NA | NA | NA | 4.00 | NA |
| South Asia | Mungbean | NA | NA | 4.50 | 3.50 | 3.50 | 3.00 |
| South Asia | Pigeonpea | NA | NA | NA | NA | 3.50 | NA |
| South Asia | Zebu cattle | 9.00 | 8.00 | NA | NA | 8.00 | 6.50 |
| South Asia | Water buffalo | 6.00 | 4.50 | NA | NA | 4.50 | NA |
| East Asia | Broomcorn Millet | 10.00 | 8.00 | NA | NA | 8.00 | NA |
| East Asia | Foxtail millet | 11.50 | 7.50 | NA | NA | 7.50 | NA |
| East Asia | Rice | 10.00 | 8.00 | 8.00 | 7.50 | 7.50 | 5.00 |
| East Asia | Soybean | 8.50 | 5.50 | NA | NA | 5.50 | 4.00 |
| East Asia | Ramie | NA | NA | NA | NA | 5.25 | NA |
| East Asia | Melon | 7.00 | 4.00 | NA | NA | 4.00 | 3.75 |
| East Asia | Pig | 12.00 | 8.50 | NA | NA | 8.50 | 6.00 |
| East Asia | Silkworm | 7.00 | 5.25 | NA | NA | 5.25 | NA |
| East Asia | Yak | NA | NA | NA | NA | 4.25 | NA |
| East Asia | Horse | 7.50 | 6.75 | 6.75 | 5.50 | 5.50 | 4.00 |
| East Asia | Bactrian Camel | NA | NA | NA | NA | 4.50 | NA |
| East Asia | Duck | 2.50 | 1.00 | NA | NA | 1.00 | NA |
| East Asia | Chicken | 6.00 | 4.00 | NA | NA | 4.00 | NA |
| New Guinea | Banana | 10.00 | 7.00 | 7.00 | 4.00 | 4.00 | NA |
| New Guinea | Taro | 10.00 | 7.00 | 7.00 | 4.00 | NA | NA |
| New Guinea | Yam | 10.00 | 7.00 | 7.00 | 4.00 | NA | NA |
| Africa and Arabia | Date palm | 7.00 | 6.00 | NA | NA | 5.00 | NA |
| Africa and Arabia | Sorghum | 8.00 | 4.00 | NA | NA | 4.00 | NA |
| Africa and Arabia | Pearl millet | NA | NA | NA | NA | 4.50 | 3.50 |
| Africa and Arabia | Fonio | NA | NA | NA | NA | 2.50 | NA |
| Africa and Arabia | Cowpea | NA | NA | NA | NA | 3.75 | NA |
| Africa and Arabia | Hyacinth bean | NA | NA | NA | NA | 3.75 | NA |
| Africa and Arabia | Rice | 3.50 | 2.00 | NA | NA | 2.00 | NA |
| Africa and Arabia | Oil palm | 9.25 | 3.50 | NA | NA | 3.50 | NA |
| Africa and Arabia | Cattle | NA | NA | 9.00 | 7.75 | 7.75 | 6.50 |
| Africa and Arabia | Donkey | 9.00 | 5.50 | NA | NA | 5.50 | 3.50 |
| Africa and Arabia | Dromedary camel | 6.50 | 3.00 | NA | NA | 3.00 | NA |
| Africa and Arabia | Guinea fowl | NA | NA | 2.50 | 1.50 | 1.50 | NA |
| North America | Squash | 6.50 | 5.00 | NA | NA | 5.00 | NA |
| North America | Sunflower | 6.00 | 4.75 | NA | NA | 4.00 | NA |
| North America | Sumpweed | 6.00 | 4.50 | NA | NA | 4.00 | NA |
| North America | Pitseed goosefoot | 4.75 | 3.75 | NA | NA | 3.75 | NA |
| Meso-america | Squash (pepo) | NA | NA | NA | NA | 10.00 | 9.50 |
| Meso-america | Maize | 10.00 | 9.00 | NA | NA | 9.00 | NA |
| Meso-america | Foxtail millet-grass | NA | NA | NA | NA | 6.00 | 4.00 |
| Meso-america | Common bean | NA | NA | NA | NA | 3.00 | NA |
| Meso-america | Avocado | NA | NA | NA | NA | 3.00 | NA |
| Meso-america | Chile pepper | NA | NA | NA | NA | 3.00 | NA |
| Meso-america | Turkey | NA | NA | NA | NA | 2.00 | NA |
| South America | Chili pepper | NA | NA | NA | NA | 6.00 | NA |
| South America | Peanut | NA | NA | NA | NA | 5.00 | NA |
| South America | Cotton | NA | NA | NA | NA | 6.00 | NA |
| South America | Coca | NA | NA | NA | NA | 8.00 | NA |
| South America | Now-minor root crops (arrowroot, leren) | NA | NA | NA | NA | 9.00 | NA |
| South America | Squash | NA | NA | NA | NA | 10.00 | NA |
| South America | Common bean | NA | NA | NA | NA | 5.00 | NA |
| South America | Lima bean | NA | NA | 8.25 | NA | 6.00 | NA |
| South America | Monioc | NA | NA | NA | NA | 7.00 | NA |
| South America | Sweet potato | NA | NA | NA | NA | 5.00 | NA |
| South America | White potato | 7.00 | 4.50 | NA | NA | 4.00 | NA |
| South America | Quinoa | 5.00 | NA | NA | NA | 3.50 | NA |
| South America | Yam | NA | NA | NA | NA | 5.50 | NA |
| South America | Llama | 10.00 | 6.00 | NA | NA | 6.00 | NA |
| South America | Alpaca | 10.00 | 5.00 | NA | NA | 5.00 | NA |
| South America | Guinea pig | NA | NA | NA | NA | 5.00 | NA |
| South America | Muscovy Duck | NA | NA | NA | NA | 4.00 | NA |
par(mar=c(5,4,6,1))
dates <- unlist(domestication_times[3:8])
hist(dates, breaks = 22, xlim=c(15,0), xlab="K years ago", col=adjustcolor("cornflowerblue", alpha= 0.5), border=adjustcolor("cornflowerblue", alpha= 0.9), main="All dates in dataset" )
mtext("This tells us about how evenly our evidence is distributed in time", 3, line=1)
hist(dates, breaks = 22, xlim=c(15,0), xlab="Thousand years ago", col=adjustcolor("cornflowerblue", alpha= 0.5), border=adjustcolor("cornflowerblue", alpha= 0.9), main="All dates in dataset with Larson(2014) date windows")
Start_of_early_window <- 12
End_of_early_window_start_of_late_window <- 8.2
End_of_late_window <- 4.2
polygon(x=c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window), y=c(0, 30, 30, 0), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x=c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window), y=c(0, 30, 30, 0), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
hist(dates, breaks = 22, xlim=c(15,0), xlab="K years ago", col=adjustcolor("cornflowerblue", alpha= 0.2), border=adjustcolor("cornflowerblue", alpha= 0.9), main="", add=TRUE)
mtext("Early Holocene", 3, line = -1, adj=.3)
mtext("Middle Holocene", 3, line= -1, adj=.6)
par(mfrow=c(2,3), mar=c(4,4,2,0))
specific_dates <- domestication_times[3:9]
Error in `[.data.frame`(domestication_times, 3:9) :
undefined columns selected
I’m creating new rows for this table, combining dates in different ways to make the CDFs below look more authentic. This makes it so that pre-ag always happens before post-ag. What I’ve done is given the later date to the earlier date when those dates are missing.
h <- which(is.na(domestication_times[,3]))
domestication_times <- cbind(domestication_times, rep(NA, length(domestication_times[,1])))
domestication_times[,9] <- domestication_times[,3]
domestication_times[h,9] <- domestication_times[h,7]
colnames(domestication_times)[9] <- "adopt exploitation date"
domestication_times[,10] <- domestication_times[,7]
domestication_times[which(is.na(domestication_times[,10])),10] <- 0
colnames(domestication_times)[10] <- "start of ag"
#save(domestication_times, file="~/Desktop/Human density and the origins of agriculture/Domestication timing larson 2014.Rdata")
I think these are best described by a cummulative distribution, showing how they accumulate over time.
for(i in 1:8){
type_number <- i
match <- domestication_times[ which(domestication_times$Region == levels(domestication_times$Region)[ type_number]), 9]
maxer <- max(match, na.rm=TRUE)
j <- ecdf(maxer-match)
print(levels(domestication_times$Region)[ type_number])
print(match)
print(j)
}
[1] "Africa and Arabia"
[1] 7.00 8.00 4.50 2.50 3.75 3.75 3.50 9.25 7.75 9.00 6.50 1.50
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 0.25, 1.25, ..., 6.75, 7.75
[1] "East Asia"
[1] 10.00 11.50 10.00 8.50 5.25 7.00 12.00 7.00 4.25 7.50 4.50 2.50 6.00
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 0.5, 2, ..., 7.75, 9.5
[1] "Meso-america"
[1] 10 10 6 3 3 3 2
Empirical CDF
Call: ecdf(maxer - match)
x[1:4] = 0, 4, 7, 8
[1] "New Guinea"
[1] 10 10 10
Empirical CDF
Call: ecdf(maxer - match)
x[1:1] = 0
[1] "North America"
[1] 6.50 6.00 6.00 4.75
Empirical CDF
Call: ecdf(maxer - match)
x[1:3] = 0, 0.5, 1.75
[1] "South America"
[1] 6.0 5.0 6.0 8.0 9.0 10.0 5.0 6.0 7.0 5.0 7.0 5.0 5.5 10.0 10.0 5.0 4.0
Empirical CDF
Call: ecdf(maxer - match)
x[1:8] = 0, 1, 2, ..., 5, 6
[1] "South Asia"
[1] 8.5 8.0 4.5 4.0 3.5 3.5 9.0 6.0
Empirical CDF
Call: ecdf(maxer - match)
x[1:7] = 0, 0.5, 1, ..., 5, 5.5
[1] "Southwest asia"
[1] 12.0 12.0 12.0 11.5 11.0 10.5 12.0 10.0 12.0 12.0 12.0 11.5 4.0
Empirical CDF
Call: ecdf(maxer - match)
x[1:6] = 0, 0.5, 1, ..., 2, 8
par(mfcol=c(2,5), mar=c(4,0,5,0))
plot(0,0, type="n", xaxt="n", xlab="", bty="n")
mtext("Percent of species that will eventually \n be domesticated in a region", 2, line=-5, cex=0.5)
plot(0,0, type="n", xaxt="n", xlab="", bty="n")
mtext("Percent of species that will eventually \n be domesticated in a region", 2, line=-5, cex=0.5)
for(i in 1:8){
type_number <- i
match <- domestication_times[ which(domestication_times$Region == levels(domestication_times$Region)[ type_number]), 9]
maxer <- max(match, na.rm=TRUE)
j <- ecdf(maxer-match)
#print(j)
plot(0,0, xlim=c(15,0), ylim=c(0,100), ylab="Percent of species that will eventually \n be domesticated in a region", xlab="Thousand years ago", main=levels(domestication_times$Region)[ type_number], type="n", yaxt="n")
x_seq <- rev(c(0,seq(0, maxer, length.out=100)))
y_seq <- 100 * (c(0, j(seq(0, maxer, length.out=100))))
lines(x_seq, y_seq, ylim=c(-1,1))
polygon(c(0, x_seq), c(0, y_seq), border=adjustcolor("cornflowerblue",alpha=1), col=adjustcolor("cornflowerblue", alpha=0.2))
if(i == 2 | i == 1)axis(2)
if(i == 3)mtext("Cummulative distribution function for the accumulation of domesticates", 3, line=3.8, col="cornflowerblue")
}
par(mfcol=c(2,5), mar=c(4,0,5,0))
plot(0,0, type="n", xaxt="n", xlab="", bty="n")
mtext("Percent of species that will eventually \n be domesticated in a region", 2, line=-5, cex=0.5)
plot(0,0, type="n", xaxt="n", xlab="", bty="n")
mtext("Percent of species that will eventually \n be domesticated in a region", 2, line=-5, cex=0.5)
for(i in 1:8){
type_number <- i
match <- domestication_times[ which(domestication_times$Region == levels(domestication_times$Region)[ type_number]), 9]
maxer <- max(match, na.rm=TRUE)
j <- ecdf(maxer-match)
#print(j)
plot(0,0, xlim=c(15,0), ylim=c(0,100), ylab="Percent of species that will eventually \n be domesticated in a region", xlab="Thousand years ago", main=levels(domestication_times$Region)[ type_number], type="n", yaxt="n")
x_seq <- rev(c(0,seq(0, maxer, length.out=100)))
y_seq <- 100 * (c(0, j(seq(0, maxer, length.out=100))))
lines(x_seq, y_seq, ylim=c(-1,1))
polygon(c(0, x_seq), c(0, y_seq), border=adjustcolor("cornflowerblue",alpha=1), col=adjustcolor("cornflowerblue", alpha=0.2))
abline(v= maxer - quantile(j)[2], col="limegreen", lwd=2)
if(i == 2 | i == 1)axis(2)
if(i == 2)mtext("25%", 3, line=3.5, adj=-1, col="limegreen")
if(i == 3)mtext("Cummulative distribution function for the accumulation of domesticates", 3, line=3.8, col="cornflowerblue")
if(i == 4)mtext("Choose a y to predict an x", 3, line=3.3, col="cornflowerblue")
break_one <- maxer
break_two <- maxer - quantile(j)[2]
polygon(x=c(break_two, break_two, break_one, break_one), y=c(0, 1, 1, 0), col=adjustcolor("cornflowerblue", alpha=0.2), border=adjustcolor("cornflowerblue",alpha=1))
lines(x=c(break_two, break_two), y=c(0,-1), col="cornflowerblue")
abline(h = 25, col="limegreen", lwd=2)
}
Make this a function. There is a choice of two methods here. At the end of this section we need to print the desision we’re passing to the later analyses.
origins
class : SpatialPolygonsDataFrame
features : 20
extent : -104.7263, 144.9554, -26.31802, 43.00233 (xmin, xmax, ymin, ymax)
coord. ref. : NA
variables : 8
no non-missing arguments, returning NAno non-missing arguments, returning NAno non-missing arguments, returning NAno non-missing arguments, returning NAno non-missing arguments, returning NAno non-missing arguments, returning NA
names : CONTINENT, SQMI, SQKM, OID, CONTINEN_2, SQMI_2, SQKM_2, OID_2
min values : C/S_Andes, NA, NA, NA, Europe, 3821854, 9898597, 7
max values : West Africa T, NA, NA, NA, Europe, 3821854, 9898597, 7
library(maps)
map()
map(origins, add=TRUE, fill=TRUE, col=adjustcolor("cornflowerblue", alpha=1))
database does not (uniquely) contain the field 'name'.
map()
d <- readPNG("Larson_origins.png")
rasterImage(d, -180, -90, 180, 110, interpolate=TRUE, col=d)
map(add=TRUE)
map(origins, add=TRUE, fill=TRUE, col=adjustcolor("cornflowerblue", alpha=1))
database does not (uniquely) contain the field 'name'.
This is obviously a bad projection fit right now.
#subset and reorder origins. This is currently done at the end of the plot but should be moved forward.
# Load data for population density
load("PopD_all_December.rdata")
PopD.ALL
class : RasterStack
dimensions : 288, 720, 207360, 18 (nrow, ncol, ncell, nlayers)
resolution : 0.5, 0.5 (x, y)
extent : -180, 180, -60, 84 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
names : fourK, fiveK, sixK, sevenK, eightK, nineK, tenK, elevenK, twelveK, thirteenK, fourteenK, fifteenK, sixteenK, seventeenK, eighteenK, ...
min values : 5.611358e-07, 1.067142e-06, 2.508241e-06, 6.317553e-06, 2.286934e-05, 7.631922e-05, 1.272693e-04, 2.118215e-04, 2.602175e-04, 3.226203e-04, 4.390267e-04, 5.572032e-04, 7.313966e-04, 8.286005e-04, 8.297062e-04, ...
max values : 2.051069, 2.013452, 2.142908, 1.888403, 1.863014, 1.880628, 1.650615, 1.678033, 1.697732, 1.499115, 1.517264, 1.443677, 1.464867, 1.453581, 1.436394, ...
# Extract data to a matrix
Pop <- values(PopD.ALL)
r <- raster(PopD.ALL, 1)
r
class : RasterLayer
dimensions : 288, 720, 207360 (nrow, ncol, ncell)
resolution : 0.5, 0.5 (x, y)
extent : -180, 180, -60, 84 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names : fourK
values : 5.611358e-07, 2.051069 (min, max)
We need to justify our decision to use a GAM over other models. This should include citations to back up those arguments.
We should make our decisions very transparent here. We should be able to justify our decision of 3 degrees of freedom over other possible values.
plot(global.means[[1]], col=adjustcolor("cornflowerblue", alpha=0.8), pch=names(global.means[[1]]), type="b", xlab="year", ylab="Density", xaxt="n")
axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
means_matrix <- matrix(rep(NA,19*20), 20, 19)
colnames(means_matrix) <- c("origin", rev(seq(4, 21, by=1)))
means_matrix[,1] <- names(global.means)
for(i in 1:20){
means_matrix[i,2:19] <- global.means[[i]]
}
kable(means_matrix, caption= "This is our world")
| origin | 21 | 20 | 19 | 18 | 17 | 16 | 15 | 14 | 13 | 12 | 11 | 10 | 9 | 8 | 7 | 6 | 5 | 4 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| W_African_Sav Means | 0.277236125574504 | 0.268146096510109 | 0.229965375382666 | 0.32523944033517 | 0.0993713185018438 | 0.0037800484090062 | 0.0288533426383454 | 0.0965336556256875 | 0.170551592391401 | 0.23617134733154 | 0.329881362382384 | 0.408474454716208 | 0.933326290541932 | 1.0205068503086 | 0.850021182238413 | 0.257169818734039 | 0.326479019859158 | 0.507805632758623 |
| Sudanic_Savan Means | 0.016989310920534 | 0.0165830329202393 | 0.0137788739758862 | 0.0187343740600394 | 0.0267628737555364 | 0.0504019933855703 | 0.111853742051633 | 0.171175628805818 | 0.229766197375977 | 0.290372573994489 | 0.393633312377065 | 0.465128867070118 | 0.952417421795436 | 1.01093274306512 | 0.893099872916669 | 0.355788155105045 | 0.546151812406506 | 0.75726422940137 |
| West Africa T Means | 0.0189861304288692 | 0.0211956898670637 | 0.0184494598528687 | 0.0122957090963052 | 0.0401148140595455 | 0.0642571026136034 | 0.122814423606705 | 0.18024456411236 | 0.238763360681693 | 0.296946782277892 | 0.393959129288487 | 0.475949031004617 | 0.875757914548295 | 0.90858669240854 | 0.874514137796506 | 0.581496815464146 | 0.745072302084471 | 0.956646116404008 |
| Ethipian plat Means | 0.0904589833106248 | 0.0903924937521069 | 0.0863100536268273 | 0.0893963604501343 | 0.101029422893253 | 0.112809652424169 | 0.165444985051315 | 0.237587289908825 | 0.335763703990818 | 0.422484059629967 | 0.563949159000512 | 0.624066354130941 | 0.940018099247958 | 0.954579684776469 | 0.920354009533087 | 0.563097081815597 | 0.604720525289939 | 0.917811255034155 |
| NA Means | 0.663114221139142 | 0.660555808084869 | 0.642678366385832 | 0.68037338350395 | 0.579197369561157 | 0.449722187452827 | 0.292062953469214 | 0.193223598846692 | 0.132827947024298 | 0.111861275509966 | 0.0710971234608467 | 0.130373982857816 | 0.671207444822536 | 0.786755748449708 | 0.739053831844885 | 0.832107257992727 | 0.655180067894376 | 0.51579662144668 |
| E_North_Ameri Means | 0.00893546567365586 | 0.0104915755544869 | 0.013253014074329 | 0.00542980780574677 | 0.033277010182152 | 0.0676074551240852 | 0.151565641055132 | 0.252887929725881 | 0.364503616958045 | 0.473544250047383 | 0.693964954332778 | 0.767791046626426 | 0.692619559646708 | 0.69478102682856 | 0.698626071405891 | 0.419432397571308 | 0.858741718110438 | 1.00985185846655 |
| New_Guinea Means | 0.280908407324953 | 0.278870321026453 | 0.267694222403951 | 0.281470118449041 | 0.252745290593119 | 0.214728303639231 | 0.167084548728992 | 0.115970302717808 | 0.089799780729113 | 0.080892066918588 | 0.0124431847502232 | 0.0401293843516986 | 0.566581929269924 | 0.604979974822395 | 0.641520404478516 | 0.955830671499772 | 0.708491079993817 | 0.662407420268634 |
| Mesoamerica Means | 0.0161950005425715 | 0.0206245055086652 | 0.02143899866462 | 0.00934611302070087 | 0.0544114903050211 | 0.103081131118678 | 0.194265350885966 | 0.278924138424292 | 0.363943243683155 | 0.43990897393249 | 0.568842547084558 | 0.615137398815427 | 0.832150878517869 | 0.835544619667127 | 0.803694000400885 | 0.462771111989956 | 0.868405844942686 | 1.00928777965156 |
| N_Lowland_SA Means | 0.197689193698132 | 0.194041814243325 | 0.190984940424971 | 0.196128712506526 | 0.196417663969857 | 0.187068023549807 | 0.160506607282115 | 0.134145108737672 | 0.107884268995178 | 0.0911509383234282 | 0.0205928310679456 | 0.0459520795917419 | 0.596005745518906 | 0.629886813716808 | 0.656340274475672 | 0.950411787925309 | 0.712319378859971 | 0.762450575454787 |
| NW_Lowland_SA Means | 0.302160479847383 | 0.3047338767083 | 0.30145828256188 | 0.306608856132229 | 0.294571575771776 | 0.269112077717038 | 0.24906509989485 | 0.206289289884621 | 0.156152857530682 | 0.120418637914445 | 0.0391309645428076 | 0.0308172782102177 | 0.297940313948749 | 0.399158518142026 | 0.500237965916322 | 0.946354452820633 | 0.691804496481964 | 0.603912759846072 |
| Sava_W_India Means | 0.0106243834199155 | 0.0114152911613111 | 0.0161244909782402 | 0.00391888342167919 | 0.0411288944114358 | 0.0799118639893852 | 0.143474824462497 | 0.194472001633401 | 0.25896453538531 | 0.328023269123674 | 0.38501994780977 | 0.438200162623175 | 0.701906926325678 | 0.887863530419312 | 0.840257997855973 | 0.463636950860515 | 0.824692040164097 | 1.00602620757196 |
| S_India Means | 0.0115629583651973 | 0.0187752993630369 | 0.021723311306329 | 0.0070646947614955 | 0.0407072830692962 | 0.0713069789187112 | 0.139798008959816 | 0.225664612623677 | 0.281610839357003 | 0.348147871179238 | 0.403540116832751 | 0.477872217545242 | 0.800247790956004 | 0.852607258819713 | 0.838151682591359 | 0.653173598407203 | 0.915365428351119 | 1.00598449001798 |
| Ganges_E_Indi Means | 0.0113365151896897 | 0.0126436772421878 | 0.0137067008861891 | 0.00706800434124777 | 0.0340469630303124 | 0.0618791589942872 | 0.107610274242122 | 0.154284668909943 | 0.214022368644172 | 0.261411688761718 | 0.354911484688783 | 0.405502150971218 | 0.701022087868857 | 0.856628657667716 | 0.857644756926133 | 0.634690322099047 | 0.846206513504704 | 1.00363122170002 |
| Chinese_loess Means | 0.00985165064501957 | 0.0104808324437268 | 0.00900683865719935 | 0.0116119701845256 | 0.0118779155381701 | 0.0441243012860027 | 0.0716263584823884 | 0.0861013791329461 | 0.149418291528035 | 0.211519249389285 | 0.266988439490975 | 0.311227904177839 | 0.707572308573447 | 0.849429071406241 | 0.906202082492626 | 0.923045544934632 | 0.9916760767357 | 0.973623792552289 |
| Japanese Means | 0.432966078309294 | 0.436721275258284 | 0.430357801245452 | 0.448302080145254 | 0.403558960792971 | 0.34383950864674 | 0.275786181477207 | 0.200371550616222 | 0.126217675525547 | 0.0867652904934615 | 0.0131109616455113 | 0.0730107508281249 | 0.809250181854901 | 0.903296939451686 | 0.864704553202643 | 0.981926278231659 | 0.863987172963444 | 0.7582788509337 |
| Lower-MiddleY Means | 0.0120866286052244 | 0.0157503645001354 | 0.0141141386633795 | 0.0117842718312604 | 0.0320108118220579 | 0.0507118754130521 | 0.091891834618743 | 0.124873052749255 | 0.158608104722965 | 0.204520350655691 | 0.243958207153347 | 0.28747052316418 | 0.680749440029052 | 0.773753745502804 | 0.84134458277978 | 0.698341733661832 | 0.907809652199586 | 1.00418953333218 |
| South trop ch Means | 0.016713449702567 | 0.0174753363663615 | 0.0151393509017418 | 0.0180283696050808 | 0.0156823173722971 | 0.0122750768667892 | 0.038639651118046 | 0.0916115579684524 | 0.114645285272606 | 0.223588081432699 | 0.297587932561125 | 0.375458659406677 | 0.805249095005994 | 0.93801021743498 | 0.945522393929682 | 0.910830218399952 | 0.966004005587459 | 0.988016065147811 |
| NA Means | 0.0506689441611688 | 0.0523003760584072 | 0.0529223473123491 | 0.0583571740189899 | 0.0488909376957515 | 0.0393647803163041 | 0.0390561263869707 | 0.133860007449605 | 0.184609429162159 | 0.197183843402241 | 0.21980312514192 | 0.355570001144045 | 0.765746252662412 | 0.921681593973495 | 0.908058913268584 | 0.886924001187376 | 0.948972971310446 | 0.951733541434419 |
| Southwes amaz Means | 0.197814786548895 | 0.196624956520459 | 0.189937961661877 | 0.190641832484773 | 0.201114079962783 | 0.192405606587582 | 0.169432074839703 | 0.14150349548783 | 0.116741484382453 | 0.108871940612822 | 0.0637616038743224 | 0.0788896738782188 | 0.603605273568561 | 0.617468071130463 | 0.651007381887985 | 0.928816410800911 | 0.927512492742812 | 0.921262845255316 |
| C/S_Andes Means | 0.0401384562944843 | 0.0397694823632042 | 0.0403745098771443 | 0.0400938475174884 | 0.0598582525765649 | 0.0714790397599819 | 0.10113089957835 | 0.115327513375983 | 0.128946585768714 | 0.144650495307302 | 0.144317872008496 | 0.186436295057576 | 0.67781043506519 | 0.704637618405442 | 0.731079856196258 | 0.895704921957057 | 0.942947005554748 | 0.961094004592374 |
#global.SD
SD_matrix <- matrix(rep(NA,19*20), 20, 19)
colnames(SD_matrix) <- c("origin", rev(seq(4, 21, by=1)))
SD_matrix[,1] <- names(global.SD)
for(i in 1:20){
SD_matrix[i,2:19] <- global.SD[[i]]
}
kable(SD_matrix, caption= "SD values")
| origin | 21 | 20 | 19 | 18 | 17 | 16 | 15 | 14 | 13 | 12 | 11 | 10 | 9 | 8 | 7 | 6 | 5 | 4 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| W_African_Sav SD | 0.172699007009174 | 0.166224154719586 | 0.138678306162886 | 0.185441293815945 | 0.0518023278157856 | 0.0176743529543277 | 0.0232721845857313 | 0.0333644331899999 | 0.0354168207585815 | 0.035199893254038 | 0.0485347261747143 | 0.0608429760322292 | 0.0242495467702773 | 0.0125228565331159 | 0.030853798105179 | 0.0816359078903485 | 0.13119632975143 | 0.136708070596833 |
| Sudanic_Savan SD | 0.0159698059987182 | 0.0157324841477218 | 0.0130990900131222 | 0.0166003275140283 | 0.015966067525986 | 0.0256027817811989 | 0.0452389607615637 | 0.0568279371893465 | 0.0652666423235462 | 0.0671560859495923 | 0.052642065744974 | 0.0523704150626269 | 0.0124936002052098 | 0.0216830502557889 | 0.0357007985386511 | 0.142302424115745 | 0.125375351923483 | 0.137098220468632 |
| West Africa T SD | 0.0167321856989961 | 0.0157899403246079 | 0.0148137837237379 | 0.0113613562356161 | 0.0196995774415615 | 0.0275689805142787 | 0.0419436547568745 | 0.050492520331069 | 0.059679399772876 | 0.0658658041127462 | 0.0828859886252774 | 0.0699049776104795 | 0.0818358041657027 | 0.092465678209583 | 0.0535103005034764 | 0.128075013406526 | 0.0795441005120006 | 0.0675621686827204 |
| Ethipian plat SD | 0.162912188924422 | 0.160411599553525 | 0.151742223557912 | 0.170525258224472 | 0.120309000695041 | 0.0832108748773215 | 0.101796496236904 | 0.150037458872254 | 0.202589035265331 | 0.229238291028049 | 0.269645889743349 | 0.253345298475953 | 0.0447447233108725 | 0.0525085647718391 | 0.0424886599682493 | 0.18245585531685 | 0.0798520618580059 | 0.102304970623811 |
| NA SD | 0.317243930545136 | 0.315343575298247 | 0.308652081354488 | 0.324003088184895 | 0.293515373848188 | 0.258502627900925 | 0.203354497097551 | 0.152404334346264 | 0.095190692866469 | 0.0644775468495577 | 0.0863609463670797 | 0.129211809891193 | 0.141308178808749 | 0.184807201581753 | 0.142524730394035 | 0.141543601124512 | 0.0976259662285693 | 0.116035132973936 |
| E_North_Ameri SD | 0.00814519801285417 | 0.00851539781161795 | 0.00929401001777065 | 0.00447318912463679 | 0.0115087080492953 | 0.0156252653953385 | 0.0225849146016072 | 0.0251984491931828 | 0.0244453007992024 | 0.0227279927606796 | 0.0320698144632027 | 0.040761082463001 | 0.0524817338555138 | 0.052946139308667 | 0.0395648579830907 | 0.0370549127832752 | 0.0262840341483891 | 0.00134154275005444 |
| New_Guinea SD | 0.0821958325250514 | 0.0828638794671788 | 0.0781541364668149 | 0.0828088029896117 | 0.078094797782274 | 0.0717879349217425 | 0.053393900283989 | 0.0494630568718273 | 0.0385500079009193 | 0.0274721883505787 | 0.0115300949179493 | 0.00800295127290903 | 0.0872756613070994 | 0.0895534405758082 | 0.0562194157260793 | 0.00552327232079952 | 0.0333048005670632 | 0.0446051787672239 |
| Mesoamerica SD | 0.0211537849845949 | 0.0229393635287474 | 0.0166418607418158 | 0.0118194786555225 | 0.0225239907442187 | 0.0352868854772045 | 0.0613620740588326 | 0.0855395320159072 | 0.103613885886215 | 0.11422095117608 | 0.129846631792337 | 0.128105208673015 | 0.0452010183656475 | 0.0556118672284279 | 0.0482593402432222 | 0.113311245060815 | 0.0396131063123546 | 0.00328698466541702 |
| N_Lowland_SA SD | 0.103306550386419 | 0.100347039338042 | 0.101714165882939 | 0.106238235672369 | 0.0977489513449625 | 0.0897789632675473 | 0.0651858939863286 | 0.0489499923305428 | 0.0381362652127659 | 0.0288166238340547 | 0.0261316061651377 | 0.0217881798432685 | 0.11032994628824 | 0.120447165832298 | 0.1068666416441 | 0.017049009533096 | 0.082446683077705 | 0.12104668727093 |
| NW_Lowland_SA SD | 0.130739696844979 | 0.132907733073528 | 0.132009693999081 | 0.133691638250333 | 0.127524836612193 | 0.113005200692562 | 0.0698516459686315 | 0.0575681624179661 | 0.059602563376245 | 0.0631030747496126 | 0.0727156805650935 | 0.0667517958417757 | 0.0781785070725703 | 0.0931558488306065 | 0.0586358262852361 | 0.00580470161060384 | 0.0452537544732729 | 0.0521503662828824 |
| Sava_W_India SD | 0.00578695471476928 | 0.00522073368988782 | 0.00613849030781516 | 0.00183399992079732 | 0.00837019252888177 | 0.0148090191021166 | 0.023478907405565 | 0.0292024405234909 | 0.033611576878901 | 0.0497449326880249 | 0.0581864938492896 | 0.0755295896901178 | 0.0829532745787329 | 0.0496220092145389 | 0.0429972870472159 | 0.0519140268417683 | 0.0445401020980192 | 0.00772933089067814 |
| S_India SD | 0.0117624884650963 | 0.0157871049212158 | 0.0118145281046149 | 0.00655509962741526 | 0.012978839356834 | 0.0126344004077163 | 0.0179081104817061 | 0.0487275050572536 | 0.0429954709946546 | 0.0282450142482082 | 0.0357676443766461 | 0.0567307976706172 | 0.038340653395058 | 0.0341238206328938 | 0.0313877845976713 | 0.0261806771552853 | 0.0276222239818159 | 0.00112405830031077 |
| Ganges_E_Indi SD | 0.00999610376800777 | 0.0101615092874938 | 0.0104876775583078 | 0.00641609933281792 | 0.0147546484547988 | 0.0234146327223517 | 0.0376291526801385 | 0.0513518755513335 | 0.0641617622544108 | 0.0700848276080226 | 0.0999052869333427 | 0.108208542272367 | 0.14721968212923 | 0.125484765455796 | 0.102598416103331 | 0.0922877014003416 | 0.0938015475022164 | 0.00253370306633808 |
| Chinese_loess SD | 0.0104576172281714 | 0.00980110612844382 | 0.00783333956532036 | 0.00756518364091116 | 0.00944388611008931 | 0.0226519110393966 | 0.0218576927676012 | 0.0148955150204466 | 0.0392830305528864 | 0.0317202170873381 | 0.046021275348067 | 0.0443694154073961 | 0.103392912258369 | 0.0993849425873263 | 0.0574571628700919 | 0.0498239200365443 | 0.00632011800334687 | 0.0230005033057195 |
| Japanese SD | 0.173874867501537 | 0.163342130431977 | 0.149019749611258 | 0.150323444589142 | 0.145901292029872 | 0.127512093574691 | 0.0700152264283969 | 0.0551989152928249 | 0.05204487446896 | 0.0482288801166012 | 0.0390583134719698 | 0.0584408591210569 | 0.0944238954809108 | 0.0443487097006404 | 0.0399553629459968 | 0.0045816516907887 | 0.0343952741808961 | 0.0528178820446734 |
| Lower-MiddleY SD | 0.0121575718392651 | 0.0135189595661818 | 0.0124669761541667 | 0.0115122244783575 | 0.0141743122719921 | 0.016363557279366 | 0.0214366417551572 | 0.0324190064000509 | 0.0392120895886192 | 0.043682275223655 | 0.0564093787608781 | 0.0562184106071083 | 0.0508853757706017 | 0.0462179408583881 | 0.0293617386404968 | 0.0480271930716818 | 0.0408102405864344 | 0.0022127237876935 |
| South trop ch SD | 0.0125153746424966 | 0.0138055463439825 | 0.0121127130053405 | 0.01426390624128 | 0.0121833901400692 | 0.010938220481384 | 0.0271378184286839 | 0.0421478968600715 | 0.0493247074658243 | 0.068195642879872 | 0.0725679264428561 | 0.0915320496460424 | 0.107401688641242 | 0.094183651074615 | 0.0669832486146481 | 0.0647028879611437 | 0.0230594488830411 | 0.0221371524026078 |
| NA SD | 0.0653732907402712 | 0.0627918419611784 | 0.0584406887359542 | 0.0651529688452993 | 0.0483153457267308 | 0.0306513656864185 | 0.0293300217951196 | 0.0650620201262001 | 0.0878013105267807 | 0.100169276861391 | 0.129586238997577 | 0.135990304726632 | 0.0574536894989798 | 0.0485457552815111 | 0.0436063435578501 | 0.0935423241238113 | 0.0428321434521962 | 0.0735379607174006 |
| Southwes amaz SD | 0.109469697830144 | 0.112448446924183 | 0.11233937395215 | 0.114499793136335 | 0.0960984647019482 | 0.0791816380668241 | 0.0489776548427169 | 0.0451477303475569 | 0.0515466212624549 | 0.0645415033910373 | 0.0879754490546639 | 0.0853931959342598 | 0.0461741195930737 | 0.0496361553947693 | 0.0443919623194883 | 0.0770832795279368 | 0.047411891055166 | 0.0660149541648784 |
| C/S_Andes SD | 0.0669388782190815 | 0.0620953064144171 | 0.0643020118362701 | 0.0647332641994508 | 0.0606662239855845 | 0.0512285580525077 | 0.0562492133276492 | 0.0626717399773927 | 0.0642146029586656 | 0.0669676599338186 | 0.084046861778514 | 0.0857024389759996 | 0.103821838349181 | 0.113118000822383 | 0.0901022957978736 | 0.0547403879595825 | 0.0711151774461137 | 0.0737932864833159 |
par(mfrow=c(5,4), mar=c(0,0,0,0))
for(j in 1:20){
plot(means_matrix[1,2:19], col=adjustcolor("cornflowerblue", alpha=0.8), pch=names(means_matrix[1,2:19]), type="n", xlab="year", ylab="Density", xaxt="n")
lines(means_matrix[j,2:19], col=adjustcolor("cornflowerblue", alpha=0.8), pch=names(means_matrix[j,2:19]), type="b", xlab="year", ylab="Density", xaxt="n")
}
axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
# Get the predctions from Population_trend script
load("prediction.RData")
# Read the polygons
origins <- readShapePoly('Origins_updated.shp')
# Extract data
cells <- do.call(rbind, sapply(per.origin, subset, select = 1))
#cells
g.means <- apply(prediction[-cells, ], 2, mean, na.rm = TRUE)
g.gams <- apply(prediction[-cells, ], 2, sd, na.rm = TRUE)
g.means2 <- apply(prediction[cells, ], 2, mean, na.rm = TRUE)
g.gams2 <- apply(prediction[cells, ], 2, sd, na.rm = TRUE)
#pdf("Global_pop_trend_comparisson.pdf", width = 25, height = 20)
par(mar = c(5, 7, 7, 5))
plot(seq(0, 1, length.out = length(time)) ~ time, col = "white", main = "GLOBAL",
xlim = c(21, 4), ylab = "Population Density (standardized)",
xlab = "Thousand of years ago", cex.lab = 1, cex.main = 1, cex.axis = 1)
down <- g.means - g.gams
up <- g.means + g.gams
lines(y = down, x = time, lty = 3, col = "gray40", lwd = 3)
lines(y = up, x = time, lty = 3, col = "gray40", lwd = 3)
lines(y = g.means, x = time, lwd = 4)
lines(y = g.means2, x = time, lwd = 3, col = "red")
down2 <- g.means2 - g.gams2
up2 <- g.means2 + g.gams2
lines(y = down2, x = time, lty = 3, col = "red", lwd = 3)
lines(y = up2, x = time, lty = 3, col = "red", lwd = 3)
polygon(cbind(c(12, 8.2, 8.2, 12, 12), c(-1, -1, 2, 2, -1)),
col = rgb(0, 1, 0, alpha = .2), border = F)
polygon(cbind(c(8.2, 4.2, 4.2, 8.2, 8.2), c(-1, -1, 2, 2, -1)),
col = rgb(.28, 0, .28, alpha = .2), border = F)
#dev.off()
# Load patricks productivity PCA data
load('Productivity_ALL.RDATA')
# Load origin shapefiles
origins <- readShapePoly('Origins_updated.shp')
origin.time.region <- c(2, 2, 1, 1, 1, 2, 2, 1, 2, 2,
2, 2, 1, 2, 2, 2, 2, 2, 2, 2) # 1 = early; 2 = middle
# Extract the data
prod.origin <- extract(Productivity.ALL, origins)
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
# Mean and SD per region
means <- lapply(prod.origin, colMeans, na.rm = TRUE)
sds <- lapply(prod.origin, sd, na.rm = TRUE)
names(means) <- origins@data$CONTINENT
ymax <- max(unlist(means))
ymin <- min(unlist(means))
time <- 4:21
# Plot
#pdf("productivity.pdf", 20, 30)
par(mfrow = c(5, 4), mar = c(2, 2, 2, 0))
for (i in 1:length(means)) {
plot(y = means[[i]], x = time, xlim = c(21, 4), ylim = c(ymin, ymax),
main = names(means)[i], cex.main = 1, cex.lab = 1, cex.axis = 1,
ylab = "Productivity (PCA axis)", xlab = "Thousand of years ago (k)",
pch = 20, lwd = 1, type = "l",
col = c("purple", "green")[origin.time.region[i]])
up <- sds[[i]] + means[[i]]
down <- means[[i]] - sds[[i]]
lines(up ~ time, lty = 2)
lines(down ~ time, lty = 2)
}
#dev.off()
a <- layout(matrix(c(
1, 1, 1, 1, 1, 1, 1, 1,
3, 6, 7, 8, 9, 10, 11, 4,
3, 5, 5, 5, 5, 5, 5, 4,
3, 12, 13, 14, 15, 16, 17, 4,
2, 2, 2, 2, 2, 2, 2, 2
), 5, 8, byrow=TRUE), width=c(1, 1, 1, 1, 1, 1, 1, 1), height=c(0.5, 1, 1.5, 1, 0.5))
layout.show(a)
frameplot <- function(){
plot(21:0,rep(0, 22), xlim=c(17,4), ylim=c(0, 2.25), type="n", xaxt="n", yaxt="n", xlab="", ylab="")
}
frameplot_bottom <- function(){
plot(21:0,rep(0, 22), xlim=c(17,4), ylim=c(-0.25, 2), type="n", xaxt="n", yaxt="n", xlab="", ylab="")
}
frameplot()
frameplot_bottom()
###################
type_number <- 8
complex_figure <- function(type_number, i, means, sds){
if(i < 6) polygon(x=c(12,12,8.2,8.2), y=c(-1,3,3,-1), col=adjustcolor("cornflowerblue", alpha=0.4), border=NA)
if(i > 5) polygon(x=c(8.2,8.2,4.2,4.2), y=c(-1,3,3,-1), col=adjustcolor("limegreen", alpha=0.4), border=NA)
match <- domestication_times[ which(domestication_times$Region == levels(domestication_times$Region)[ type_number]), 9]
maxer <- max(match, na.rm=TRUE)
j <- ecdf(maxer-match)
print(j)
x_seq <- rev(c(0,seq(0, maxer, length.out=100)))
y_seq <- -c(0, j(seq(0, maxer, length.out=100)))
#lines(x_seq, y_seq, type="l", ylim=c(-1,1))
#polygon(c(0, x_seq), c(0, y_seq), border="black", col=adjustcolor("cornflowerblue", alpha=0.5))
#abline(v= maxer - quantile(j)[2])
break_one_1 <- maxer
break_two_1 <- maxer - quantile(j)[2]
# polygon(x=c(break_two_1, break_two_1, break_one_1, break_one_1), y=c(0, 1, 1, 0), col=adjustcolor("cornflowerblue", alpha=0.5), border=NA)
match <- domestication_times[ which(domestication_times$Region == levels(domestication_times$Region)[ type_number]), 10]
maxer <- max(match, na.rm=TRUE)
j <- ecdf(maxer-match)
print(j)
x_seq <- rev(c(0,seq(0, maxer, length.out=100)))
y_seq <- 2+c(0, j(seq(0, maxer, length.out=100)))
#lines(x_seq, y_seq)
#polygon(c(0, x_seq), c(2, y_seq), border="black", col=adjustcolor("limegreen", alpha=0.5))
break_one_2 <- maxer
break_two_2 <- maxer - quantile(j)[2]
# polygon(x=c(break_two_2, break_two_2, break_one_2, break_one_2), y=c(1, 2, 2, 1), col=adjustcolor("limegreen", alpha=0.5), border=NA)
#abline(v=11)
type <- 1
if(type == 1){
x <- c(means[[i]] , means[[i]] + abs(sds[[i]]), means[[i]] - abs(sds[[i]]))
scaled <- scale(x , center=FALSE)
meanss <- scaled[1:18]
sdss_plus <- scaled[19:36]
sdss_minus <- scaled[37:54]
#abline(v=10, col="red")
length(scaled)
#lines(4:21, means[[i]] + sds[[i]])
#polygon(x=c(4:21, 21:4), y=c(sdss_plus, rev(sdss_minus)), col=adjustcolor("firebrick", alpha=1), border="white")
polygon(x=c(21:4,4:21), y=c(sdss_plus, rev(sdss_minus)), col=adjustcolor("firebrick", alpha=1), border="white")
}
if(type == 2){
x <- c(means[[i]] , means[[i]] + abs(sds[[i]]), means[[i]] - abs(sds[[i]]))
scaled <- x + 1 #scale(x , center=FALSE)
meanss <- scaled[1:18]
sdss_plus <- scaled[19:36]
sdss_minus <- scaled[37:54]
#abline(v=10, col="red")
length(scaled)
#lines(4:21, means[[i]] + sds[[i]])
polygon(x=c(4:21, 21:4), y=c(sdss_plus, rev(sdss_minus)), col=adjustcolor("firebrick", alpha=1), border="white")
}
if(type == 3){
x <- c(means[[i]] , means[[i]] + abs(sds[[i]]), means[[i]] - abs(sds[[i]]))
scaled <- x #scale(x , center=FALSE)
meanss <- scaled[1:18]
sdss_plus <- scaled[19:36]
sdss_minus <- scaled[37:54]
#abline(v=10, col="red")
length(scaled)
#lines(4:21, means[[i]] + sds[[i]])
polygon(x=c(4:21, 21:4), y=c(sdss_plus, rev(sdss_minus)), col=adjustcolor("firebrick", alpha=1), border="white")
}
means_long_y <- c(1,1,1,1,1, meanss)
means_long_x <- c(0:4, 4:21)
break_one <- break_one_2
break_two <- break_two_2
# polygon(x=c(break_one, break_one, 22, 22), y=c(1, 2, 2, 1), col=adjustcolor("white", alpha=0.8), border=NA)
# polygon(x=c(break_two, break_two, break_one, break_one), y=c(1, 2, 2, 1), col=adjustcolor("white", alpha=0), border=NA)
# polygon(x=c(-1,-1, break_two , break_two), y=c(1.9, 3.1, 3.1, 1.9), col=adjustcolor("white", alpha=0.8), border=NA)
#abline(v= break_one, col="white")
#abline(v= break_two, col="white")
break_one <- break_one_1
break_two <- break_two_1
# polygon(x=c(break_one, break_one, 22, 22), y=c(0, 1, 1, 0), col=adjustcolor("white", alpha=0.8), border=NA)
# polygon(x=c(break_two, break_two, break_one, break_one), y=c(0, 1, 1, 0), col=adjustcolor("white", alpha=0), border=NA)
# polygon(x=c(-1,-1, break_two , break_two), y=c(-1.1, .1, .1, -1.1), col=adjustcolor("white", alpha=0.8), border=NA)
#abline(v= break_one, col="white")
#abline(v= break_two, col="white")
#lines(x=c(break_one_2, break_one_2), y=c(1,3), col="white")
#lines(x=c(break_one_1, break_one_1), y=c(1,-1), col="white")
#lines(x=c(break_two_2, break_two_2), y=c(1,3), col="white")
#lines(x=c(break_two_1, break_two_1), y=c(1,-1), col="white")
#lines(4:21, meanss)
lines(21:4, meanss)
}
frameplot()
complex_figure(7, 1, global.means, global.SD)
Empirical CDF
Call: ecdf(maxer - match)
x[1:7] = 0, 0.5, 1, ..., 5, 5.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:4] = 0, 3.5, 4, 4.5
axis(1)
axis(2)
quartz(width=8, height=8)
layout(matrix(c(
1, 1, 1, 1, 1, 1, 1, 1,
3, 6, 7, 8, 9, 10, 11, 4,
3, 5, 5, 5, 5, 5, 5, 4,
3, 12, 13, 14, 15, 16, 17, 4,
2, 2, 2, 2, 2, 2, 2, 2
), 5, 8, byrow=TRUE), width=c(1, 1, 1, 1, 1, 1, 1, 1), height=c(0.5, 1, 1.5, 1, 0.5))
par(mar=c(0,0,0,0))
# 1-4 label margins
blankplot <- function(){
plot(0,0, xlim=c(4,21), ylim=c(1, 1.25), bty="n", type="n", xaxt="n", yaxt="n", xlab="", ylab="")
}
blankplot()
blankplot()
blankplot()
blankplot()
origins <- readShapePoly('Origins_updated.shp')
origin.time.region <- c(2, 2, 1, 1, 1, 2, 2, 1, 2, 2,
2, 2, 1, 2, 2, 2, 2, 2, 2, 2) # 1 = early; 2 = middle
as.character(origins$CONTINENT)
[1] "W_African_Sav" "Sudanic_Savan" "West Africa T" "Ethipian plat" "Fertile_Cresc"
[6] "E_North_Ameri" "New_Guinea" "Mesoamerica" "N_Lowland_SA" "NW_Lowland_SA"
[11] "Sava_W_India" "S_India" "Ganges_E_Indi" "Chinese_loess" "Japanese"
[16] "Lower-MiddleY" "South trop ch" "Chinese_loess" "Southwes amaz" "C/S_Andes"
#subset_order <- c(1, 2, 3, 5, 6, 8, 9, 10, 11, 12, 17, 18)
subset_order <- c(8, 10, 9, 5, 18, 7, 6, 20, 1, 2, 13, 16)
origins_subset <- origins[subset_order,]
origins_subset$CONTINENT
[1] Mesoamerica NW_Lowland_SA N_Lowland_SA Fertile_Cresc Chinese_loess New_Guinea
[7] E_North_Ameri C/S_Andes W_African_Sav Sudanic_Savan Ganges_E_Indi Lower-MiddleY
19 Levels: C/S_Andes Chinese_loess E_North_Ameri Ethipian plat ... West Africa T
d <- readPNG("earth.png")
png(file=paste("40962.png",sep=""),width=2000,height=1000, bg="transparent")
par(mar=c(0,0,0,0))
plot(seq(-180, 180, length.out = 19), seq(-90, 90, length.out = 19), type="n",xlim=c(-180, 180),ylim=c(-90, 90), xaxt="n")
rasterImage(d, -180, -90, 180, 90, interpolate=TRUE, col=d)
polygon(x=c(-180,-180, 180,180), y=c(-90, 90, 90, -90), col=adjustcolor("white", alpha=0.1))
#rasterImage(d, -13.5, -13.5, 375, 375, interpolate=TRUE, col=d)
plot(origins_subset, add=TRUE, col=adjustcolor("white", alpha=.8), xaxt="n", border="white", lwd=4) #still need to reproject!!!
dev.off()
quartz
2
d <- readPNG("40962.png")
dim(d)
[1] 1000 2000 4
par(mar=c(0,0,0,0))
plot(0:360,0:360,type="n",xlim=c(20,360),ylim=c(65,295), yaxt="n", xaxt="n")
rasterImage(d, -28.5, -13.5, 388, 375, interpolate=TRUE, col=d)
axis(2, label=seq(-90, 90, length.out = 19), at=seq(1, 360, length.out = 19), las=1)
mtext("latitude", 2, line=4, at=180)
abline(h=seq(1, 360, length.out = 19), col=adjustcolor("grey10", alpha= 0.4), lwd=1)
abline(h=180, col=adjustcolor("white", alpha= .5), lwd=1)
load('PopD_all_December.rdata')
# Extract the data
prod.origin <- extract(PopD.ALL, origins_subset)
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
library(matrixStats)
# Mean and SD per region
means <- lapply(prod.origin, colMeans, na.rm = TRUE)
sds <- lapply(prod.origin, colSds, na.rm = TRUE)
## new values from Bruno's GAM model (produced in script called Population_Trend_per_y.R)
means <- global.means
sds <- global.gams
names(means) <- origins_subset@data$CONTINENT
ymax <- max(unlist(means))
ymin <- min(unlist(means))
time <- 4:21
#plot(origins)
#means[[1]] +
#sds[[1]]
#scale(as.numeric(means[[1]]), center=FALSE)
name_vector <- as.character(origins_subset@data$CONTINENT)
for(i in 1:12){
if(i > 6){frameplot()}else{frameplot_bottom()}
## customize polygons for each graph
if(i == 1){ #mesoamerica #values from Larson
complex_figure(3, i, means, sds)
}
#########
if(i == 2 ){ #NW lowlands SA #values from Larson
complex_figure(6, i, means, sds)
}
#########
if( i == 3){ #NW lowlands SA #values from Larson
complex_figure(6, i, means, sds)
}
#########
if(i == 4){ #Fertile crescent aka Southwest asia #values from Larson
complex_figure(8, i, means, sds)
}
#########
if(i == 5){ #loess plateau #values from Larson
complex_figure(2, i, means, sds)
}
#########
if(i == 6){ #new guinea #values from Larson
complex_figure(4, i, means, sds)
}
#########
if(i == 7){ #Eastern N.A. #values from Larson
complex_figure(5, i, means, sds)
}
#########
if(i == 8){ #Andes #values from Larson
complex_figure(6, i, means, sds)
}
#########
if(i == 9){ #W. African Sav #values from Larson
complex_figure(1, i, means, sds)
}
#########
if(i == 10){ #Sudanic sav #values from Larson
complex_figure(1, i, means, sds)
}
#########
if(i == 11){ #Ganges #values from Larson
complex_figure(7, i, means, sds)
}
#########
if(i == 12){ #loess #values from Larson
complex_figure(2, i, means, sds)
}
#lines(4:21, means[[i]])
abline(h = 1, col=adjustcolor("forestgreen", alpha=.5), lty=2)
# add axes to some locations
if(i == 1 | i == 7){axis(2, at=seq(0,2, by=0.25), label=seq(0,2, by=0.25), las=1)}
if(i == 6 | i == 12){axis(4, at=seq(0,2, by=0.25), label=seq(0,2, by=0.25), las=1)}
#if(i == 6 | i == 12){axis(4, at=seq(2,3, by=0.25), label=seq(0,1, by=0.25), las=1)
# axis(4, at=seq(-1,0, by=0.25), label=rev(seq(0,1, by=0.25)), las=1)
# }
if(i > 6){axis(1)} else{axis(3)}
# add text
if(i < 7){polygon(x=c(-30, -30, 30, 30), y=c(-0.1, -0.5, -0.5, -0.1), col="black")
mtext(name_vector[i], 1, line=-1.2, col="white", cex=0.5)}
if(i > 6){polygon(x=c(-30, -30, 30, 30), y=c(2.1, 2.5, 2.5, 2.1), col="black")
mtext(name_vector[i], 3, line=-1.2, col="white", cex=0.5)}
# add axis labels
if(i == 1 | i == 7){mtext("scaled density potential", 2, line=4, at=1)}
if(i == 3){mtext("Thousand years before present", 3, line=3.5, at =5)}
if(i == 9){mtext("Thousand years before present", 1, line=3.5, at =5)
}
}
Empirical CDF
Call: ecdf(maxer - match)
x[1:4] = 0, 4, 7, 8
Empirical CDF
Call: ecdf(maxer - match)
x[1:5] = 0, 1, 4, 7, 8
Empirical CDF
Call: ecdf(maxer - match)
x[1:8] = 0, 1, 2, ..., 5, 6
Empirical CDF
Call: ecdf(maxer - match)
x[1:9] = 0, 1, 2, ..., 6, 6.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:8] = 0, 1, 2, ..., 5, 6
Empirical CDF
Call: ecdf(maxer - match)
x[1:9] = 0, 1, 2, ..., 6, 6.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:6] = 0, 0.5, 1, ..., 2, 8
Empirical CDF
Call: ecdf(maxer - match)
x[1:8] = 0, 0.5, 0.75, ..., 5, 7
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 0.5, 2, ..., 7.75, 9.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:9] = 0, 0.5, 1, ..., 4.5, 7.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:1] = 0
Empirical CDF
Call: ecdf(maxer - match)
x[1:2] = 0, 4
Empirical CDF
Call: ecdf(maxer - match)
x[1:3] = 0, 0.5, 1.75
Empirical CDF
Call: ecdf(maxer - match)
x[1:3] = 0, 1, 1.25
Empirical CDF
Call: ecdf(maxer - match)
x[1:8] = 0, 1, 2, ..., 5, 6
Empirical CDF
Call: ecdf(maxer - match)
x[1:9] = 0, 1, 2, ..., 6, 6.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 0.25, 1.25, ..., 6.75, 7.75
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 2.25, 2.75, ..., 5.75, 6.25
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 0.25, 1.25, ..., 6.75, 7.75
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 2.25, 2.75, ..., 5.75, 6.25
Empirical CDF
Call: ecdf(maxer - match)
x[1:7] = 0, 0.5, 1, ..., 5, 5.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:4] = 0, 3.5, 4, 4.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 0.5, 2, ..., 7.75, 9.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:9] = 0, 0.5, 1, ..., 4.5, 7.5
saveToPDF <- function(...) {
d = dev.copy(pdf,...)
dev.off(d)
}
saveToPNG <- function(...) {
d = dev.copy(png,...)
dev.off(d)
}
## Try them out
saveToPDF("my.pdf", height=8,width=8)
quartz
2
saveToPNG("my.png", height=8, width=8, units="in", res=300)
quartz
2
dev.off()
null device
1